14 Area Data V

Area Data V

NOTE: You can download the source files for this book from here. The source files are in the format of R Notebooks. Notebooks are pretty neat, because the allow you execute code within the notebook, so that you can work interactively with the notes.

In the previous chapter, you learned about how to use local spatial statistics for exploratory spatial data analysis.

If you wish to work interactively with this chapter you will need the following:

• An R markdown notebook version of this document (the source file).

• A package called geog4ga3.

Learning Objectives

In this practice, you will:

1. Practice how to estimate regression models in R.
2. Learn about autocorrelation as a model diagnostic.
3. Learn about variable transformations.
4. Use autocorrelation analysis to improve regression models.

Suggested Readings

• Bailey TC and Gatrell AC [-@Bailey1995] Interactive Spatial Data Analysis, Chapter 7. Longman: Essex.
• Bivand RS, Pebesma E, and Gomez-Rubio V [-@Bivand2008] Applied Spatial Data Analysis with R, Chapter 9. Springer: New York.
• Brunsdon C and Comber L [-@Brunsdon2015R] An Introduction to R for Spatial Analysis and Mapping, Chapter 7. Sage: Los Angeles.
• O’Sullivan D and Unwin D [-@Osullivan2010] Geographic Information Analysis, 2nd Edition, Chapter 7. John Wiley & Sons: New Jersey.

Preliminaries

As usual, it is good practice to clear the working space to make sure that you do not have extraneous items there when you begin your work. The command in R to clear the workspace is rm (for “remove”), followed by a list of items to be removed. To clear the workspace from all objects, do the following:

rm(list = ls())

Note that ls() lists all objects currently on the worspace.

Load the libraries you will use in this activity:

library(tidyverse)
#library(rgdal)
#library(geosphere)
library(ggmap)
library(geosphere)
#library(gridExtra)
library(sf)
library(plotly)
library(spdep)
#library(crosstalk)
library(geog4ga3)

Next, read an object of class sf (simple feature) with the census tracts of Hamilton CMA and some selected population variables from the 2011 Census of Canada:

data(Hamilton_CT)

The sf object can be converted into a SpatialPolygonsDataFrame object for use with the spdedp package:

Hamilton_CT.sp <- as(Hamilton_CT, "Spatial")

Regression Analysis in R

This section provides a refresher on linear regression, before reviewing the estimation of regression models in R.

A linear regression model posits relationships between an outcome, called a dependent variable, and one or more covariates, called independent variables. Regression models capture statistical relationships, not causal relationships, but often the causality is implied by the choice of independent variables.

This is the form of a linear regression model: \[ y_i = \beta_0 + \sum_{j=1}^k{\beta_jx_{ij}} + \epsilon_i \] where \(y_i\) is the dependent variable and \(x_ij\) (\(j=1,...,k\)) are the independent variables. The coefficients \(\beta\) are not known, but can be estimated from the data. And \(\epsilon_i\) is a term called a residual (or error), because it is the difference between the systematic term of the model and the value of \(y_i\): \[ \epsilon_i = y_i - \bigg(\beta_0 + \sum_{j=1}^k{\beta_jx_{ij}}\bigg) \]

Estimation of a linear regression model is the procedure used to obtain values for the coefficients. This typically involves defining a loss function that needs to be minimized. In the case of linear regression, a widely used estimation procedure is least squares. This procedure allows a modeler to find the coefficients that minimize the sum of squared residuals. In very simple terms, the protocol is as follows: \[ \text{Find the values of }\beta\text{ that minimize }\sum_{i=1}^n{\epsilon_i^2} \]

For this procedure to be valid, there are a few assumptions that need to be satisfied, including:

1. The functional form of the model is correct.

2. The independent variables are not collinear; this is often diagnosed by calculating the correlations among the independent variables, with values greater than 0.8 often being problematic.

3. The residuals have a mean of zero: \[ E[\epsilon_i|X]=0 \]

4. The residuals have constant variance: \[ Var[\epsilon_i|X] = \sigma^2 \text{ }\forall i \]

5. The residuals are independent, that is, they are not correlated among them: \[ E[\epsilon_i\epsilon_j|X] = 0 \text{ }\forall i\neq j \]

The last three assumptions ensure that the residuals are random. Violation of these assumptions is often a consequence of a failure in the first two (i.e., the model was not properly specified, and/or the residuals are not exogenous).

When all these assumptions are met, the coefficiens are said to be BLUE: Best Linear Unbiased Estimates - a desirable property because we wish to be able to quantify the relationships between covariates without bias.

The basic command for multivariate linear regression in R is lm(), for “linear model”. This is the help file of this function:

?lm #Remember that we can search the definition of a function by means of placing a question mark in front of the function itself. 

Lets see how to estimate a model using this function.

To illustrate this, we will use the example of population density gradients. These gradients are representations of the variation of population density in cities, and are of interest because they are related to land rent, urban form, and commuting patterns, among other things (see accompanying reading for more information).

Urban economic theory suggests that population density declines with distance from the central business district of a city, or its CBD. This leads to the following model, where the population density at \(i\) is a function of the distance of \(i\) to the CBD. Since this is likely a stochastic process, we allow for some randomness by means of the residuals: \[ P_i = f(D_i) + \epsilon_i \]

Lets add distance to the CBD as a covariate to your dataframe. We will use Jackson Square as the CBD of Hamilton:

xy_cbd <- c(-79.8708, 43.2584)
?c

Convert these coordinates to sf, with EPSG code 4326 for longitude and latitude:

#xy_cbd <- st_as_sf(xy_cbd, coords = c("lon", "lat"), crs = 4326)

xy_ct <- coordinates(spTransform(Hamilton_CT.sp, CRSobj = "+proj=longlat +datum=WGS84 +no_defs"))
#we use 'spTransform' as a method for map projection.

?coordinates

Given these coordinates, the function sf::st_distance can be used to calculate the great circle distance between the centroids of the census tracts and Hamilton’s CBD. Call this dist2cbd.sl, i.e., straight line distance to CBD in a straight:

dist2cbd.sl <- distGeo(xy_ct, xy_cbd)
?distGeo

Add the distance to Hamilton_CT for analysis:

Hamilton_CT$dist.sl <- dist2cbd.sl

Regression analysis is implemented in R by means of the lm function. The arguments of the model include an object of type “formula” and a dataframe. Other arguments include conditions for subsetting the data, using sampling weights, and so on.

A formula is written in the form y ~ x1 + x2, and more complex expressions are possible too, as we will see below. For the time being, the formula is simply POP_DENSIT ~ dist.sl:

model1 <- lm(formula = POP_DENSITY ~ dist.sl, data = Hamilton_CT)
summary(model1)

Call:
lm(formula = POP_DENSITY ~ dist.sl, data = Hamilton_CT)

Residuals:
    Min      1Q  Median      3Q     Max 
-3841.2 -1338.3  -177.1   950.9 10009.1 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 4405.15414  250.16355  17.609  < 2e-16 ***
dist.sl       -0.17989    0.02418  -7.439 3.63e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1892 on 186 degrees of freedom
Multiple R-squared:  0.2293,    Adjusted R-squared:  0.2251 
F-statistic: 55.33 on 1 and 186 DF,  p-value: 3.631e-12

The value of the function is an object of class lm that contains the results of the estimation, including the coefficients with their significance diagnostics, and the coefficient of multiple determination, among other items.

Notice how the coefficient for distance is negative (and significant). This indicates that population density declines with increasing distance: \[ P_i = f(D_i) + \epsilon_i = 4405.15414 - 0.17989D_i + \epsilon_i \]

Autocorrelation as a Model Diagnostic

Let quickly explore the fit of the model. The scatterplot is of the actual population density and the distance to CBD, whereas the blue line is the regression model:

ggplot(data = Hamilton_CT, aes(x = dist.sl, y = POP_DENSITY)) + 
  geom_point() +
  geom_abline(slope = model1$coefficients[2], 
              intercept = model1$coefficients[1], 
              color = "blue", size = 1) +
  geom_vline(xintercept = 0) + geom_hline(yintercept = 0)

Clearly, there is a fair amount of noise (the scatter of the dots around the regression line). In this case, the regression line captures the general trend of the data, but underestimates most of the high population density areas, and overestimates most of the low population areas.

If the pattern of under- and over-estimation is random (i.e., the residuals are random), that would indicate that the model successfully retrieved all the systematic pattern. If the pattern is not random, there is a violation of assumption of independence.

Lets add the residuals of the model to your dataframes:

Hamilton_CT$model1.e <- model1$residuals

And now we can create a map of the residuals, with red indicating negative residuals and blue positive:

  plot_ly(Hamilton_CT) %>% #Recall 'plot.ly' is an interactive mapping tool
    add_sf(color = ~(model1.e > 0), colors = c("red", "dodgerblue4")) #You can add characteristics to a plotly map by means of the 'add_sf' function. This example adds colours to represent positive and negative residuals

Does the distribution of residuals look random?

Fortunately, we already have the tools to help us with this question. Lets create a set of spatial weights:

Hamilton_CT.w <- nb2listw(poly2nb(Hamilton_CT.sp)) #Recall that we can create a lsit of spatial weights of nearest neighbours using 'nb2listw' 
?poly2nb

With this, we can calculate Moran’s I:

moran.test(Hamilton_CT$model1.e, Hamilton_CT.w)

    Moran I test under randomisation

data:  Hamilton_CT$model1.e  
weights: Hamilton_CT.w    

Moran I statistic standard deviate = 9.7447, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.395382878      -0.005347594       0.001691092 

This supports the visual inspection of the map, since the test of hypothesis rejects the null hypothesis of independence at a very high level of confidence (note the extremely small value of p).

Spatial autocorrelation, as mentioned above, is a violation of a key assumption of linear regression, and likely the consequence of a model that was not correctly specified, either because the functional form was incorrect (e.g., the relationship was not linear), or there are missing covariates.

Lets explore the first of these possibilities by means of variable transformations.

Variable Transformations

The term linear regression refers to the linearity in the coefficients. Variable transformations allow you to consider non-linear relationships between covariates, while still preserving the linearity of the coefficients.

For instance, a possible transformation of the variable distance could be as inverse distance: \[ f(D_i) = \beta_0 + \beta_1\frac{1}{D_i} \]

Lets create a new covariate that is the inverse distance:

Hamilton_CT <- mutate(Hamilton_CT, invdist.sl = 1/dist.sl)
#Recall that mutate adds to existing variables while preserving those that already exist. We are mutating Hamilton census trat information to create inverse distance

Use inverse distance as the covariate in the model:

model2 <- lm(formula = POP_DENSITY ~ invdist.sl, data = Hamilton_CT)
summary(model2)

Call:
lm(formula = POP_DENSITY ~ invdist.sl, data = Hamilton_CT)

Residuals:
   Min     1Q Median     3Q    Max 
 -6763  -1375    -52   1108   9675 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)    2299.6      164.4  13.988  < 2e-16 ***
invdist.sl  2260521.5   342039.7   6.609 3.97e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1940 on 186 degrees of freedom
Multiple R-squared:  0.1902,    Adjusted R-squared:  0.1858 
F-statistic: 43.68 on 1 and 186 DF,  p-value: 3.965e-10

As the scatterplot below shows (the model is the blue line), we can capture a non-linear relationship, which does a better job of describing the high density of tracts close to the CBD, but unfortunately is a poor description of density almost everywhere else:

ggplot(data = Hamilton_CT, aes(x = dist.sl, y = POP_DENSITY)) + 
  geom_point() +
  stat_function(fun=function(x)2299.6 + 2260521.5/x, geom="line", color = "blue", size = 1) +
  geom_vline(xintercept = 0) + geom_hline(yintercept = 0)

Add the residuals of this model to the dataframe for further examination:

Hamilton_CT$model2.e <- model2$residuals

If we calculate Moran’s I, you will notice that the coefficient is lower than for the previous model but still significant, wich means that the residuals are not random:

moran.test(Hamilton_CT$model2.e, Hamilton_CT.w)

    Moran I test under randomisation

data:  Hamilton_CT$model2.e  
weights: Hamilton_CT.w    

Moran I statistic standard deviate = 8.8236, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.358132814      -0.005347594       0.001696970 

The literature on population density gradients suggests other non-linear transformations, including: \[ f(D_i) = exp(\beta_0)exp(\beta_1x_i) \]

This function is no longer linear in the coefficients (since all the coefficientsare transformed by the exponential). Fortunately, there is a simple way of changing this to a linear expression, by taking the logarithm on both sides of the equation: \[ ln(P_i) = \beta_0 + \beta_1x_i + \epsilon_i \]

Lets create these covariates:

Hamilton_CT <- mutate(Hamilton_CT, lnPOP_DEN = log(POP_DENSITY))

And estimate the model:

model3 <- lm(formula = lnPOP_DEN ~ dist.sl, data = Hamilton_CT)
summary(model3)

Call:
lm(formula = lnPOP_DEN ~ dist.sl, data = Hamilton_CT)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.5857 -0.3395  0.2970  0.6897  2.4224 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  8.465e+00  1.588e-01  53.294  < 2e-16 ***
dist.sl     -1.161e-04  1.536e-05  -7.561 1.77e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.202 on 186 degrees of freedom
Multiple R-squared:  0.2351,    Adjusted R-squared:  0.231 
F-statistic: 57.18 on 1 and 186 DF,  p-value: 1.77e-12

ggplot(data = Hamilton_CT, aes(x = dist.sl, y = POP_DENSITY)) + 
  geom_point() +
  stat_function(fun=function(x)exp(8.465 - 0.0001161 * x), 
                geom="line", color = "blue", size = 1) +
  geom_vline(xintercept = 0) + geom_hline(yintercept = 0)

As before, add the residuals of the model to the dataframe for further examination:

Hamilton_CT$model3.e <- model3$residuals

While this latest model provides a somewhat better fit, there is still systematic under- and over-prediction, as seen in the map below (red are negative residuals and blue are positive):

  plot_ly(Hamilton_CT) %>%
    add_sf(color = ~(model3.e > 0), colors = c("red", "dodgerblue4"))

Moran’s I as well strongly suggests that the residuals are not independent:

moran.test(Hamilton_CT$model3.e, Hamilton_CT.w)

    Moran I test under randomisation

data:  Hamilton_CT$model3.e  
weights: Hamilton_CT.w    

Moran I statistic standard deviate = 8.0156, p-value = 5.478e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.325928113      -0.005347594       0.001708056 

A Note about Spatial Autocorrelation in Regression Analysis

Spatial autocorrelation was originally seen as a problem in regression analysis. It is not difficult to see why, after testing three models in this chapter.

I prefer to view spatial autocorrelation as an opportunity for discovery. For instance, the models above all seem to struggle capturing the large variations in population density between the central parts of the city and the suburbs of Hamilton. Perhaps this could be due to a regime change, or in other words, the presence of an underlying process that operates somewhat differently in different parts parts of the city. The last model, for instance, suggests that the close proximity of Burlington might have an effect.

The analysis that follows is somewhat more advanced, but serves to illustrate the idea of spatial autocorrelation as a tool for discovery.

Lets begin by creating local Moran maps to identify potential spatial regimes:

localmoran.map(Hamilton_CT, Hamilton_CT.w, "POP_DENSITY", by = "TRACT")

Examination of the map above, suggests that there are possibly three regimes: a CBD (“HH” and significant tracts), Suburbs (“LL” and significant tracts), and Other (not significant tracts). Based on this, we will create two indicator variables, one for census tracts in the CBD and another for census tracts in the Suburbs. An indicator variable takes values of 1 or zero, depending on whether a condition is true. For instance, all census tracts in the CBD will take the value of 1 in the CBD indicator variable, and all others will be zero.

Begin by computing the local statistics:

POP_DEN.lm <- localmoran(Hamilton_CT$POP_DENSITY, listw = Hamilton_CT.w)

Next, identify the type of tract based on the spatial relationships (i.e., “HH”, “LL”, or “HL/LH”). After that, identify as CBD all tracts for which Type is “HH” and the p-value is less than or equal to 0.05. Likewise, identify as Suburb all tracts for which Type is “LL” and the p-value is also less than or equal to 0.05:

df_msc <- transmute(Hamilton_CT,
                    TRACT = TRACT,
                    Z = (POP_DENSITY - mean(POP_DENSITY)) / var(POP_DENSITY),
                    SMA = lag.listw(Hamilton_CT.w, Z),
                    Type = factor(ifelse(Z < 0 & SMA < 0, "LL",
                                         ifelse(Z > 0 & SMA > 0, "HH", "HL/LH"))))
df_msc <- cbind(df_msc, POP_DEN.lm)

CBD <- ifelse(df_msc$Type == "HH" & df_msc$`Pr.z...0` < 0.05, 1, 0)
Suburb <- ifelse(df_msc$Type == "LL" & df_msc$`Pr.z...0` < 0.05, 1, 0)

Add the indicator variables to the dataframe:

Hamilton_CT$CBD <- CBD
Hamilton_CT$Suburb <- Suburb

The model that I propose to estimate is a variation of the last non-linear specification, but with regime breaks: \[ ln(P_i) = \beta_0 + \beta_1x_i + \beta_2CBD_i + \beta_3Suburb_i + \beta_4CBD_ix_i + \beta_5Suburb_ix_i + \epsilon_i \]

Since the indicator variables for CBD and Suburb take values of zero and one, effectively we have the following: \[ ln(P_i)=\Bigg\{\begin{array}{l l} (\beta_0 + \beta_2) + (\beta_1 + \beta_2)x_i + \epsilon_i \text{ if census tract } i \text{ is in the CBD}\\ (\beta_0 + \beta_3) + (\beta_1 + \beta_5)x_i + \epsilon_i \text{ if census tract } i \text{ is in the Suburbs}\\ \beta_0 + \beta_1x_i + \epsilon_i \text{ otherwise}\\ \end{array} \]

Notice how the model now allows for different slopes and intercepts for observations in different parts of the city. Estimate the model:

model4 <- lm(formula = lnPOP_DEN ~ dist.sl + CBD * dist.sl + Suburb * dist.sl,
             data = Hamilton_CT)
summary(model4)

Call:
lm(formula = lnPOP_DEN ~ dist.sl + CBD * dist.sl + Suburb * dist.sl, 
    data = Hamilton_CT)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.4523 -0.2128  0.1806  0.4907  2.3234 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)     7.943e+00  1.571e-01  50.577  < 2e-16 ***
dist.sl        -3.248e-05  1.561e-05  -2.081  0.03880 *  
CBD             9.536e-01  3.788e-01   2.518  0.01267 *  
Suburb         -9.535e-01  5.089e-01  -1.874  0.06260 .  
dist.sl:CBD    -2.730e-05  1.410e-04  -0.194  0.84669    
dist.sl:Suburb -1.006e-04  3.556e-05  -2.828  0.00521 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9384 on 182 degrees of freedom
Multiple R-squared:  0.5435,    Adjusted R-squared:  0.5309 
F-statistic: 43.33 on 5 and 182 DF,  p-value: < 2.2e-16

The model is a much better fit (see the the coefficient of multiple determination).

Lets examine the residuals:

Hamilton_CT$model4.e <- model4$residuals
plot_ly(Hamilton_CT) %>%
  add_sf(color = ~(model4.e > 0), colors = c("red", "dodgerblue4"))

It is not clear from the visual inspection that the residuals are independent, but this can be tested as usual by means of Moran’s I coefficient:

moran.test(Hamilton_CT$model4.e, Hamilton_CT.w)

    Moran I test under randomisation

data:  Hamilton_CT$model4.e  
weights: Hamilton_CT.w    

Moran I statistic standard deviate = 0.74498, p-value = 0.2281
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.024432414      -0.005347594       0.001597957 

Based on the results, we fail to reject the null hypothesis, and can be confident that the residuals are likely random.

The following figure illustrates the final model:

fun.1 <- function(x)exp(7.943 + 0.9536 - (0.00003248  + 0.0000273) * x) #CBD
fun.2 <- function(x)exp(7.943 - 0.9535 - (0.00003248 + 0.0001006) * x) #Suburb
fun.3 <- function(x)exp(7.943 - 0.00003248 * x) #Other

ggplot(data = Hamilton_CT, aes(x = dist.sl, y = POP_DENSITY)) +
  geom_point() +
  geom_point(data = filter(Hamilton_CT, CBD == 1), color = "Red") +
  geom_point(data = filter(Hamilton_CT, Suburb == 1), color = "Blue") +
  stat_function(fun= fun.1, 
                geom="line", size = 1, aes(color = "CBD")) +
  stat_function(fun=fun.2, 
                geom="line", size = 1, aes(color = "Suburb")) +
  stat_function(fun=fun.3, 
                geom="line", size = 1, aes(color = "Other")) +
  scale_color_manual(values = c("red", "black", "blue"), 
                     labels = c("CBD", "Other", "Suburb")) +
  geom_vline(xintercept = 0) + geom_hline(yintercept = 0)